%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              Finate Sample performance                                    %
%                    Quantile Regression vs Kernel Nonparametric                            %
%                    Curse of Dimensionality and Misspecification                           %
%                                    Experiment 4                                           %                                   
%                         True model: V = g_0 + 1*x1 + g_2*x2                               %
%                           Estimated model: V = g_0 + g_2*x2                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
                                    
clear all

L=100; % sample size
N=2; % number of bidders
lgrid=200; % grid for quantile levels
alpha=linspace(0,1,lgrid); % quantile levels
m=size(alpha,2);
gamma0=@(alpha)(-0.1*(log(1- (alpha/exp(1))))); % intercept (function of quantile level)
gamma2=@(alpha)(1-exp(-alpha));
psi=N*(alpha.^(N-1))-(N-1)*(alpha.^N); % Psi function given N=2

toler=0.01; % parameter of tolerance for estimation
nsim=1000; % number of simulations

% Auxiliary for summulative sum
MSE_res_l = 0;
MSE_res_np_l = 0;
MSE_ER_l = 0;
MSE_ER_np_l = 0;

% Simulation Exercise
for j=1:nsim

    % Data Generating Process   
     x=[ones(L,1),lognrnd(0,0.5,L,1),exprnd(1,L,1)];
     xmed=median(x);
     xmiss=x(:,[1,3]);
     xmissmed=median(xmiss);
     p=size(xmiss,2);
    
     for k=1:m
        V(:,k)=gamma0(alpha(k)) + x(:,2) + gamma2(alpha(k))*x(:,3); % True Private Value Quantile Function
     end


     % Private Values are random draws from V with replacement and winning bid is WB=V(2:I)
     for l=1:L
        v_d=randsample(V(l,:),N,true);
        v_dd=sort(v_d, 'descend');
     % Winning bids
        wb(l,:)=v_dd(:,2);
     end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    % Quantile Regression Estimation %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for k=1:m

        gammai=[gamma0(psi(k)),gamma2(psi(k))]; % initial guess for gamma
        [~,~,gammas]=MM(ones(L,1)*psi(k), p, wb, xmiss, toler, gammai'); % gamma_hat
        V_adjk=xmiss*gammas; % Adjusted V
        ER1(:,k)=V_adjk*(N*(N-1)*(alpha(k)^(N-2))*(1-alpha(k))); % estimate of the integrand in eq (2.9) for each k
        ER2(:,k)=V_adjk*N*(alpha(k)^(N-1))*(1-alpha(k)); % estimate of the first two terms in eq (2.9) for each k
        ER1_true(:,k)=V(:,k)*(N*(N-1)*(alpha(k)^(N-2))*(1-alpha(k))); % true value of the integrand in eq (2.9) for each k
        ER2_true(:,k)=V(:,k)*N*(alpha(k)^(N-1))*(1-alpha(k)); % true value of the first two terms in eq (2.9) for each k

    end
     
     % Computation of optimal expected payoff and reservation price for each auction l:
    for l=1:L
        
        % Computation of expected payoff when screening level is in [Alpha(1),Alpha(m-1)]
        for k=1:m-1   
            intmed=trapz(alpha(k:m),ER1(l,k:m),2); % integrate ER1 from Alpha(k) to Alpha(m)
            ER_QR(:,k)=ER2(l,k)+intmed; % estimate of the expected payoff
            intmed_true=trapz(alpha(k:m),ER1_true(l,k:m),2); % integrate ER1_true from Alpha(k) to Alpha(m)
            ER_QR_true(:,k)=ER2_true(l,k)+intmed_true; % true value of the expected payoff
        end
        
        % Computation of expected payoff when screening level is Alpha(m)
        ER_QRm=[ER_QR,ER2(l,m)]; % estimate
        ER_QRm_true=[ER_QR_true,ER2_true(l,m)]; % true value
        
        % Maximization of expected payoff (estimate)
        ER_QRmax(l,:)=max(ER_QRm,[],2); % max of expected payoff
        [~,indx]=ismember(ER_QRmax(l,:),ER_QRm);
        qindx=indx; % index in Alpha that maximize estimated expected payoff
        q_rstar=alpha(qindx); % optimal screening level
        psi_q_rstar=N*(q_rstar^(N-1))-(N-1)*(q_rstar^N); % psi function used in the computation of the optimal reservation price
        gammai=[gamma0(psi_q_rstar),gamma2(psi_q_rstar)]; % initial guess for estimation
        [~,~,g_q_rstar]=MM(ones(L,1)*psi_q_rstar, p, wb, xmiss, toler, gammai'); % gamma_hat corresponding to the optimal screening level
        % Optimal Res Price
        V_rstar(l,:)=xmiss(l,:)*g_q_rstar; % optimal reservation price
        
        % Maximization of expected payoff (true)
        ER_QRmax_true(l,:)=max(ER_QRm_true,[],2); % max of true expected payoff
        [~,indx_true]=ismember(ER_QRmax_true(l,:),ER_QRm_true);
        qindx_true=indx_true; % index in Alpha that maximize true expected payoff
        V_rstar_true(l,:)=V(l,qindx_true); % optimal reservation price under the true value of V
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% 

        %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        % Nonparametric Estimation %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%

        s=sqrt(diag(cov(x(:,3)))); % standard deviation of x
        swb = sqrt(cov(wb)); % standard deviation of wb
        % Grid for wb 
        mawb=max(wb)+0.2;
        WB=linspace(0,mawb,m);
        WB=WB';
        % Generate bandwidth sequence
        h = [swb,s]*2.34*(L^(-1/5));
        hv = [h(:,1)/3,2*h(:,1)/3,h(:,1),4*h(:,1)/3,3*h(:,1);h(:,2)/3,2*h(:,2)/3,h(:,2),4*h(:,2)/3,3*h(:,2)]'; % choices of bandwidth

        % Nonparametric estimation of the winning bid density and distribution function and private value distribution 
        for b=1:length(hv)
        for k=1:m
            Fwb(l,k)=cdf_psi(WB(k,:),xmiss(l,2),hv(b,:),wb,xmiss(:,2)); % winning bid distribution
            fwb(l,k)=pdf_psi(WB(k,:),xmiss(l,2),hv(b,:),wb,xmiss(:,2)); % winning bid density
            Fv(l,k) = 1 - sqrt(1 - Fwb(l,k)); % private value distribution
        end

        fexp = WB'.*fwb(l,:); % nonparametric estimation of the integrand in the seller expected payoff 


        for i=1:length(WB)-1
            cexp = trapz(WB(i:m,:),fexp(i:m),2); % integrate fexp from WB(i) to WB(m-1)
            ER_NPi(i,:) =  WB(i,:)*(Fwb(l,i) - Fv(l,i)^2) + cexp; % seller's expected payoff evaluated in the grid of WB
        end

        ER_NPm = WB(m,:)*(Fwb(l,m) - Fv(l,m)^2); % expected payoff when wb = upper boundary of the winning bids support
        ER_NP = [ER_NPi;ER_NPm]; % nonparametric estimate of the expected payoff
        
        % Maximization of expected payoff 
        ER_NP_max(l,b)=max(ER_NP); % max of expected payoff
        [~,indx]=ismember(ER_NP_max(l,b),ER_NP);
        r_star(l,b)=WB(indx); % optimal reservation price
        end

    end
    
    % Mean Square Error (quantile regression estimation)
    MSE_res_l = MSE_res_l + (V_rstar - V_rstar_true).^2; % MSE of the optimal reservation price
    MSE_ER_l = MSE_ER_l + (ER_QRmax - ER_QRmax_true).^2; % MSE of the expected payoff
    
    % Mean Square Error (nonparametric estimation)
    MSE_res_np_l = MSE_res_np_l + (r_star - (V_rstar_true*ones(1,length(hv)))).^2; % MSE of the optimal reservation price
    MSE_ER_np_l = MSE_ER_np_l + (ER_NP_max - (ER_QRmax_true*ones(1,length(hv)))).^2; % MSE of the expected payoff

    j
     

end

% Mean Squared Error of the Estimators
RMSE_res_qr = sqrt((1/(nsim*L))*sum(MSE_res_l));
RMSE_ER_qr = sqrt((1/(nsim*L))*sum(MSE_ER_l));
RMSE_res_np = sqrt((1/(nsim*L))*sum(MSE_res_np_l,1));
RMSE_ER_np = sqrt((1/(nsim*L))*sum(MSE_ER_np_l,1));

